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ABSTRACT 

Differentiation  arithmetic  is  an  ordered-pair  arithmetic  which  evaluates  both  the  value 
and  derivative  of  functions  defined  by  formulas  or  subroutines,  without  symbolics  or  ap¬ 
proximations.  As  in  the  case  of  complex  arithmetic,  multiplication  and  division  are  defined 
in  terms  of  several  real  operations.  Algorithms  are  given  for  evaluation  of  these  operations 
with  the  same  accuracy  as  real  multiplication  and  division,  that  is.  to  the  closest  floating¬ 
point  number.  The  same  kind  of  optimal  implementation  is  described  for  Taylor  arithmetic, 
which  permits  calculation  of  Taylor  coefficients  of  arbitrary  order  for  functions  defined  by 
formulas  or  subroutines.  •»  . 
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SIGNIFICANCE  AND  EXPLANATION 


It  is  commonly  believed  that  evaluation  of  derivatives  of  a  function  requires  symbolic 
manipulation  or  numerical  approximation.  However,  for  functions  defined  by  formulas  or 
computer  subroutines,  evaluation  of  the  function  using  an  ordered-pair  arithmetic  called 
differentiation  arithmetic  yields  exact  values  of  both  the  function  and  its  derivative  at  a 
given  value  of  the  variable.  In  order  to  implement  differentiation  arithmetic  with  maximum 
accuracy  on  a  computer,  special  algorithms  have  to  be  used  for  the  derivative  components 
of  multiplication  and  division,  since  these  are  defined  by  several  real  operations,  and  hence 
are  subject  to  more  roundoff  error  than  the  corresponding  real  operations.  Algorithms 
are  given  which  implement  multiplication  and  division  in  differentiation  arithmetic  to  the 
same  accuracy  as  real  arithmetic,  each  component  of  the  result  suffering  a  single  roundoff 
error.  This  allows  calculation  of  derivatives  with  the  same  level  of  roundoff  error  as  function 
values.  For  higher  derivatives  (more  precisely,  Taylor  coefficients),  an  arithmetic  of  (n-e  1)- 
tuples  called  Taylor  arithmetic  can  be  used  in  the  same  way  as  differentiation  arithmetic 
for  the  first  derivative.  The  optimal  implementation  of  Taylor  arithmetic  is  indicated, 
based  on  well-known  algorithms  in  the  theory  of  optimal  computer  arithmetic. 


with  MRC,  and  not  with  the  author  of  this  report. 
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OPTIMAL  IMPLEMENTATION  OF  DIFFERENTIATION  ARITHMETIC 


L.  B.  Rail 


1.  Automatic  Differentiation 

The  process  of  automatic  differentiation  8-  consists  essentially  of  the  evaluation  of  func¬ 
tions  defined  b\  formulas  or  finite  algorithms  using  rules  for  operations  and  representations 
of  variables  and  constants  pertaining  to  differentiation  arithmetic  ill]  rather  than  ordinary 
real  arithmetic.  This  process  is  similar  to  the  evaluation  of  a  function  /  :  R  — ♦  R  given  by 
a  formula,  such  as 


(1.1) 


/(*) 


(x  -  1)  •  (x  3) 
x  +  2 


in  complex  instead  of  real  arithmetic  to  obtain  the  complex  value  f(z)  of  the  function 
/  :  C  — »  C  defined  by  (1.1)  with  the  real  variable  x  replaced  by  the  complex  variable  z. 
This  reinterpretation  of  (1.1)  as  a  complex  function  can  be  expressed  by  saying  that  its 
symbols  have  been  “overloaded”  by  their  corresponding  meanings  as  complex  operators, 
variables  and  constants. 

Evaluation  of  (1.1)  in  differentiation  arithmetic  is  carried  out  by  a  similar  overload¬ 
ing.  Like  complex  arithmetic,  differentiation  arithmetic  is  an  ordered-pair  arithmetic,  with 
elements  of  the  form  U  =  (u,u'),  where  u,u'  €  R.  Complex  numbers,  of  course,  are  also 
usually  represented  on  a  computer  by  pairs  z  —  (x, y),  where  x,y  €  R  and  z  =  x-f-  ty.  The 
rules  of  differentiation  arithmetic  are  as  follows: 


(1.2) 


V  -  V  =  (u,u')  —  (u.t/)  =  (u  +  v,u'  +  t/), 


(1.3) 


U  -  V  =  (u ,u')  -  ( v ,  t/)  =  (u  -  v,  u'  -  v'). 
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(1.4) 


U  •  V  —  (u,u')  •  (r, vj  =  (ti  •  v,u  ■  v'  +  v  •  u;), 


U/V  -  (u,  u')/(v.  v')  -  (u/r,  (ti#  -  u  •  v '  jv)iv).  v  0. 


The  operations  inside  the  parentheses  in  the  above  definitions  are,  of  course,  real 
operations  on  real  numbers.  It  is  easy  to  recognize  that  the  first  components  of  the  results 
follow  the  rules  of  real  arithmetic,  while  the  second  components  implement  the  rules  for 
differentiation  of  the  results  of  the  corresponding  operations,  given  that  u,v  represent 
values  of  functions,  and  u',  v'  values  of  their  derivatives  at  some  real  x.  Since 

„  v  dx  dc 

(,'6)  S  =  S=0' 

where  x  denotes  the  independent  variable  and  c  a  constant,  it  follows  that  if  x  in  (1.1)  is 
replaced  by  X  =  (i,  1),  and  the  constants  1,2,3  by  (1,0),  (2,0)  and  (3,0)  respectively,  then 
evaluation  using  the  rules  of  differentiation  arithmetic  gives 


/(X)  =  /((x,l))  =  (/(x),/'(x)), 


where  /(x)  and  f(x)  denote  respectively  the  value  of  the  real  function  /  :  R  — >  R  and  its 
derivative  /'  :  R  — ♦  R  at  the  real  value  x. 

Denoting  the  set  R2  with  the  rules  of  arithmetic  ( 1 .2)— ( 1 .5)  by  D,  the  result  (1.6) 
can  be  viewed  as  overloading  the  symbols  of  (l.l)  to  obtain  the  corresponding  mapping 
/  :  D  — >  D  in  the  same  way  the  complex  extension  of  /  was  obtained.  From  an  algebraic 
standpoint,  D  is  a  division  ring  with  an  identity  element  [  1 1 J . 

This  process  is  called  automatic  differentiation  to  distinguish  it  from  symbolic  differ¬ 
entiation,  which  results  in  a  formula  for  f(x)  (which  would  then  have  to  be  evaluated),  and 
numerical  differentiation,  which  produces  only  an  approximate  value  for  f'(x).  Automatic 
differentiation  requires  only  the  formula  or  algorithm  for  /,  and  obtains  exact  values  for 
/(x)  and  f'(x).  This  combines  the  advantages  of  symbolic  and  numerical  differentiation. 


while  avoiding  their  disadvantages  (an  extra  evaluation  process  for  /'( x)  in  the  symbolic 
case,  or  approximation  error  if  the  differentiation  is  done  numerically). 

It  is  interesting  to  note  that  the  rules  for  addition  and  subtraction  in  D  are  the  same  as 
for  the  complex  numbers  C.  while  the  rules  (1.4)  and  (1.5)  for  multiplication  and  division 
in  D  are  respectively  simpler  than  the  corresponding  rules  in  C.  The  identity  element  of 
addition  in  D  is  0  -  (0.0),  while  the  identity  for  multiplication  is  1  —  (1,0).  Elements 
of  D  of  the  form  (0,  u')  are  called  nonunits.  The  nonunits  are  divisors  of  zero,  because 
(0,  u')  •  (0,  v')  —  (0,0)  for  arbitrary  u',v'.  However,  if  (u,u')  •  (v,t/)  =  (0,0)  and  u  ^  0, 
then  (u,v')  =  (0,0)  [ll]. 

Automatic  differentiation  is  not  limited  to  rational  functions.  In  the  notation  intro¬ 
duced  above,  the  chain  rule  of  calculus  can  be  expressed  as 

(18)  fW )  =  /((«, »'))  =  (/(«). «'  •  /'(«)), 

where  /  :  R  — »  R  is  a  differentiable  real  function  and  U  6  D.  Thus,  for  example,  the 
exponential  function  is  defined  by 

e"  =  e{u’u,)  =  (eu,u'-eu), 

and  so  on  [8],  [ll].  However,  the  object  of  this  paper  is  to  focus  on  optimal  computer 
implementation  of  the  four  arithmetic  operations  (l.2)-(1.5)  in  D. 


2.  Computer  Arithmetic 


The  discussion  here  follows  the  general  theory  of  computer  arithmetic  as  given  in  [6j.  On 
a  computer,  one  works  with  a  finite  subset  S  C  R.  which  is  called  the  set  of  floating-point 
numbers.  The  elements  of  S  may  vary  from  one  machine  to  another,  but  S  is  assumed 
to  satisfy  certain  natural  conditions  in  each  case,  namely  that  S  contains  0,1,  and  the 
negative  -s  of  each  s  C  S  !6  . 

In  addition  to  the  real  numbers,  most  mathematical  systems  which  appear  in  numerical 
mathematics  can  be  constructed  as  Cartesian  products  of  a  suitable  number  of  copies  of 
R,  with  appropriate  definitions  of  operations  in  terms  of  real  arithmetic.  For  example, 
the  complex  numbers  C  and  elements  of  D  can  be  considered  to  be  ordered  pairs  of 
real  numbers,  that  is,  elements  of  R  x  R.  The  usual  representation  of  these  systems  on  a 
computer  is  by  corresponding  Cartesian  products  of  S.  Thus,  ordered  pairs  of  floating-point 
numbers  (elements  of  S  x  S)  are  the  elements  of  the  computer  representations  CS  and 
DS  of  C  and  D,  respectively.  Representations  MS  C  M  of  other  numerical  mathematical 
systems  M  in  terms  of  floating-point  numbers  can  be  constructed  in  a  similar  way.  It  is 
assumed  that  MS  contains  the  negatives  of  each  of  its  elements,  and  the  identity  elements 
for  addition  and  multiplication  of  elements  of  M. 

In  mathematical  systems  M  which  are  constructed  by  taking  Cartesian  products  of 
R,  each  element  m  €  M  consists  of  components  mi,  each  of  which  is  a  real  number.  Thus, 
it  is  possible  to  introduce  a  partial  ordering  <  of  M  componentwise,  that  is,  for  m,n  €  M, 

(2.1)  m  <  n  <=>  m,  <  n,  Vi. 


The  same  definition  provides  a  partial  ordering  of  the  floating-point  representation  MS  of 
M  because  MS  C  M.  This  allows  the  definition  of  rounding  from  M  to  MS.  which  is  a 


mapping  □:  M  -♦  MS  with  the  properties  of  a  semimorphism  6): 


(2  2) 


□  m  =  m,  Vm  6  MS, 

min  jm  <D  n,  V  comparable  m, n  f  M, 

l  (  m)  ---  3]  m,  Vm  £  M. 


In  the  case  of  a  rounding  from  R  to  S  satisfying  (2.2),  the  second  condition  implies 
that  there  are  no  floating  point  numbers  (elements  of  S)  between  r  and  □  r  for  each  r  £  R. 
In  this  sense,  the  rounding  is  said  to  be  of  maximal  accuracy.  The  case  most  commonly 
considered  of  rounding  from  R  to  S  is  the  rounding  of  an  element  r  £  R  to  the  closest 
element  of  S  to  r,  with  ties  broken  by  a  rule  which  satisfies  the  third  condition  of  (2.2). 
This  rounding,  which  will  be  denoted  by  Oi  *s  said  to  be  of  maximum  accuracy.  The  ideas 
of  maximal  and  maximum  accuracy  of  roundings  apply  componentwise  to  mathematical 
systems  M  which  are  products  of  R.  Thus,  there  will  be  no  elements  of  S  between  m, 
and  □  m,  in  the  first  case,  and  Om»  is  the  closest  element  of  S  to  each  component  m,  of 
m  €  M  in  the  second. 

The  concept  (2.2)  of  a  semimorphism  provides  a  rational  way  in  which  to  define 
operations  in  MS  on  the  basis  of  the  corresponding  operations  in  M.  Suppose,  for  example, 
that  *  denotes  a  binary  operator  in  M,  for  example,  one  of  the  arithmetic  operators 
+  i .  Then,  the  corresponding  operator  in  MS,  denoted  by  E,  is  defined  by 


(2.3)  m[±\n  =□  (m  ★  n),  Vm,n€  MS. 

Implementation  of  the  operator  ★  in  MS  in  this  way  is  said  to  be  optimal  with  respect  to 
the  rounding  □. 


3.  Implementation  of  Differentiation  Arithmetic 

A  method  for  the  optimal  implementation  of  differentiation  arithmetic  with  respect  to 
the  rounding  Q  :  D  — ■  DS  of  maximum  accuracy  will  now  be  described.  First  of  all, 
suppose  that  the  rounding  :  R  — >  S  and  corresponding  optimal  arithmetic  operators 
are  available.  Then,  operators  can  be  defined  for  U,V  <E  DS  by 

U  \'  :=-  (u  ©  v,  u'  0  t/), 

U  -  l  :=  (u  Q  v,  u'  ©  v'), 


(3.1) 


j  U  ■  V  :=  (u  Q  v,  u  ©  v'  ©  v  ©  u'). 


[  U/V  :=  (u  0  v,  (u1  QuQv' 0v)  0v),  v  £  0. 

The  above  is  an  example  of  the  “vertical”  definition  of  computer  arithmetic  [6]  in 
a  mathematical  system  in  which  arithmetic  is  defined  componentwise  in  terms  of  real 
arithmetic.  Even  though  the  rounding  O  :  R  - 1 "  S  of  maximum  accuracy  is  used,  one  has 
in  general  that 


Q(u  •  v  4  v  ■  u)  ±  u  ©  v'  ©  v  ©  u', 

(3.2) 

0((u'  -  u  ■  v' /v)/v)  (u1  ©  u  ■  v'  0  v)  0  V. 

Hence,  the  computer  implementation  (3.1)  of  differentiation  arithmetic,  while  easy  to  pro¬ 
gram  for  applications  [8],  [9P  [lOj,  is  suboptimal.  Special  algorithms  are  required  for  the 
implementation  of  multiplication  and  division  in  DS. 

Optimal  implementation  of  multiplication  can  be  accomplished  by  use  of  the  optimal 
scalar  product  of  real  vectors,  which  is  required  by  the  general  theory  of  computer  arith¬ 
metic.  That  is,  if  x,y  6  Sn,  say  x  -  [x\ ,  x2, . . .  ,x„)  and  y  -  (y\,y2-,  ■  •  • ,  yn),  then  their 
scalar  product  is  the  real  number 


n 

(3.3)  i*  y  y„ 

i=i 

sometimes  also  called  the  dot  product,  or  inner  product,  of  x  and  y.  Optimal  implemen¬ 
tation  of  this  product  with  respect  to  the  rounding  requires  the  operator  (*  defined 
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(3.4) 


Vx,y  €  S". 


x©y:=  Ql  2^x«-yiJ- 

According  to  the  definition  of  .this  means  that  the  closest  element  of  S  to  the  sum 
of  products  on  the  right  sides  of  (3.3)  and  (3.4)  has  to  be  calculated.  It  turns  out  that 
implementation  of  ©  is  not  difficult  il  ,  and  this  operator  is  already  available  in  Pascal-SC 
and  the  IBM  product  ACR1TH.  Thus,  an  optimal  implementation  of  multiplication  in  DS 
with  respect  to  ©  >s  given  by 

(3.5)  U  V  :=  (uQv,  (u,v)  ©  (t/,u')),  V£/,KeDS. 

Division  in  DS  requires  a  more  careful  approach.  Setting  W  =  ( w,w ')  =  U /V ,  one 
can  write 

,  v  ■  u'  —  u  •  v1 

(3.6)  w  —  u/v,  w  - - ; - , 

as  in  elementary  calculus.  The  optimal  rounding  O w  =  u  0  v  of  the  first  component  w 
of  W  presents  no  problem.  The  numerator  of  w'  can  be  computed  to  the  closest  floating¬ 
point  number  by  using  the  optimal  scalar  product  (u,v)  ©  (-t>',  u'),  sis  in  the  case  of 
multiplication,  and  0(1’2)  =  vQv.  Hence,  the  numerator  and  denominator  of  w'  can  each 
be  computed  to  the  closest  floating-point  number,  but  their  quotient  can  differ  slightly  from 
O'  f 2).  There  is  also  a  practical  problem,  which  also  arises  in  the  algorithm  for  optimal 
complex  division  (13],  namely,  that  overflow  or  underflow  can  occur  in  the  computation  of 
the  numerator  or  denominator,  even  though  the  exact  result  w'  can  be  represented  with 
maximum  accuracy  by  an  element  of  S.  This  problem  can  be  handled  by  the  introduction  of 
appropriate  scale  factors,  followed  by  calculation  of  the  scaled  result  to  maximum  accuracy, 
and  then  examination  of  the  scale  factors  to  determine  if  overflow  or  underflow  will  actually 
occur  [  131.  If  not,  then  a  good  approximation  to  C  w'  is  obtained  in  this  way. 


In  order  to  correct  w'„  if  necessary,  the  method  of  iterative  refinement  of  the  solution 
of  a  linear  systems  of  equations  is  used  [2  ,  3:,  12  .  Note  that  W  =  U/V  is  equivalent  to 
V  •  IV  -  f\  which  gives  the  lower  triangular  system  of  equations 


(3.7) 


V  ■  W  -  ti, 

I  /  I 

v  ■  ir  -  ?•  •  w  —  u  , 


for  the  components  w,w'  of  W  in  terms  of  the  components  of  U,  V .  Thus,  given  approxi¬ 
mations  wo,Wo  to  the  solutions  of  (3.7),  for  w  =  w0  +  t  and  w’  —  w’0  -f  e\  one  has 


(  =  (u  -  v  ■  wu)/ v. 


(3.8) 

e1  =  (u1  -  v  •  w'0  -  v'  ■  Wo  -  v'  •  ()/v. 

Since  u>o  has  been  computed  to  maximum  accuracy,  c  is  negligible  with  respect  to  u)o,  but 
could  effect  the  value  of  c'.  Again,  the  numerator  of  c'  is  calculated  by  use  of  the  optimal 
scalar  product,  giving 


(3.9)  w\  :=  Wq  ®  ((u',  -v, -v1)  ®  (l,u>o,c))  0  v. 

In  this  case,  one  correction  is  sufficient  [2],  and  one  has 


(3.10) 


Qw  =  w0,  Qw' =  w\. 


In  summary,  an  algorithm  for  optimal  division  in  DS  is  as  follows: 
1  .  If  u  -  0  then 


(3.11)  w o  :=  0,  w\  u' 0  c, 
else  compute 

(3.12)  u\,:-uCi\  ((u.  v)  ®  (~  t>\  u'))  0(e  ■'  r). 


scaling  the  numerator  and  denominator  of  ir'  if  necessary  to  avoid  over  or  underflow'; 
if  u  ',  -  S.  then  calculate  tr',  by  (3.9); 


2  .  Set 


(3.13)  UjV  :=  (iwo.ii.-;)  e  DS. 

Thus,  it  is  possible  to  implement  the  operators  *  in  DS  with  maximum  accuracy 
using  the  corresponding  operators  *  in  S  and  the  optimal  scalar  product  0.  This  is 
an  example  of  the  “horizontar  definition  of  computer  arithmetic  in  the  floating-point 
system  DS  corresponding  to  the  mathematical  system  D  of  differentiation  arithmetic. 
The  suboptimal  definition  (3.1)  of  computer  arithmetic  in  DS  will  result  in  derivatives 
being  computed  with  a  greater  number  of  rounding  errors  than  function  values,  while  the 
definitions  of  addition  and  subtraction  in  (3.1),  multiplication  by  (3.5),  and  division  by 
(3.13)  will  result  in  equal  numbers  of  rounding  errors  for  both  function  and  derivative 
evaluations. 

The  techniques  of  Bohm  |2],  [3]  for  evaluation  of  arithmetic  expressions  with  maximum 
accuracy  can  also  be  extended  from  S  to  DS,  because  derivatives  of  arithmetic  operations 
are  again  defined  by  arithmetic  expressions  in  (1.2)-(1.5). 
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V  v*  v  vi^  ^  w  - - - - - 


4.  Taylor  Arithmetic 

The  definition  of  differentiation  arithmetic  given  in  §1  can  be  generalized  immediately  to 
give  values  of  higher  derivatives,  or,  more  conveniently,  Taylor  coefficients 

(4.1)  Uk  =  «(/c)(*)  •  k  -  0.1 - - 

of  the  function  u(x  -  h),  adopting  the  usual  convention  for  k  =  0  that  =  u(x)  [8j.  The 
corresponding  Taylor  polynomial  u„(x  4-  h)  of  degree  n  in  h  is  uniquely  defined  by  the 
vector 

(4.2)  U  =  (u0,  Ul,...,  tin) 

in  Rra+1.  As  in  the  case  of  D,  arithmetic  operations  are  defined  for  such  vectors  to 
obtain  a  mathematical  system  T„.  Addition  and  subtraction  are  the  corresponding  vector 
operations: 


(4.3) 


U  ±  V  =  (u0  ±  v0,u i  ±  wj, . . .  ,u„  ±  r„). 


Multiplication  and  division  are  derived  from  the  corresponding  operations  between  Taylor 
polynomials,  with  the  result  truncated  to  degree  n.  For  multiplication, 


(4.4) 


W  =  U  V,  wk  =  ^u,  k  =  0,1, - n. 


1  =  0 


In  the  case  of  division.  W  =  U jV  is  equivalent  to  V  •  IV  —  U ,  so  (4.4)  can  be  used  to  obtain 


{k  - 1  N 

Uk  w’  '  Vk~'  I  / v",  k  ~  0, 1 , .  •  • , n,  v0 -/ 0. 


i  =  0 


One  sees  immediately  that  R  -  To.  D  =  Tj.  For  1  <  n  <  +oc,  the  Tn  are  division 
rings  with  identity,  the  same  as  D.  For  n  =  oo,  the  arithmetic  defined  by  (4.3)  (4.5)  is 
simply  the  arithmetic  of  formal  power  series  (see  5],  Chapter  1),  and  T-^  is  an  integral 


domain.  It  follows  immediately  from  the  definition  of  the  arithmetic  operations  in  T„ 


that  evaluation  of  a  rational  function  u  :  R  —  R  at  X  =  (x,h,0, . . .  ,0)  €  Tn  gives 
U  =  u(X)  =  (uo,Ui , . . .  ,  u„)  £  T„,  where  the  u k  are  the  Taylor  coefficients  (4.1)  [8'. 

Automatic  generation  of  Taylor  series  is  not  limited  to  rational  functions.  As  in  the 
case  of  D.  it  is  also  possible  to  define  standard  functions  in  Tn  by  use  of  well-known 
recurrence  relations  for  their  Taylor  series  expansions  ;8i.  For  example,  for  W  -  e'  = 

exp{6'}  =  (u’t,.  «'i . u>„),  one  has 

k- l  ^  . 

(4.6)  iv ft  =  eu“,  Wk  =  Wi  Uk~"  *  =  1»2 

«=0 

The  floating-point  system  corresponding  to  Tn  will  be  denoted,  as  before,  by  TnS. 
Optimal  implementations  of  +  » — » •  with  respect  to  the  rounding  O  are  easily  performed 
using  the  corresponding  operators  ®  in  S  and  the  optimal  scalar  product  ©.  For  example, 
defining 

(4.7)  Uk  =  (u0,ui,...,u*),  V_*  =  (w*,  vjfc_j ,  —  ,  v0). 

one  has 

(4.8)  U  V  :=  {u0Qv0,Ui  ©V_„t/2®V_2 . Un  ®  V.n). 

Just  as  in  DS,  division  is  the  only  arithmetic  operator  which  requires  a  special  algo¬ 
rithm  for  its  optimal  implementation  in  T„S.  A  good,  but  suboptimal,  implementation 
can  be  made  immediately  on  the  basis  of  (4.5)  [4].  Define  the  vectors  V'_k  £  Sfc  by 

(4.9)  V'_k  =  (vfc,wfc_i,...,vi). 

Then,  (4.5)  can  be  approximated  closely  by 

(4.10)  Wi,  :=  ur ,  0  Vn,  wk  '■=  (^Uk  ~  Wk-  i  ©  Vjfj  O  t’,,.  k  -  1.2 . n. 

An  algorithm  for  optimal  division  in  T„S  can  be  developed  as  in  DS  by  noting  that 
W  -  U /V  is  equivalent  to  V  -  W  =  V .  which  in  turn  is  equivalent  to  the  lower  triangular 


system  of  equations 


V0  •  tV(i  =  Uo, 

Vj  •  w(l  -r  Vr,  ■  Wj  =  ut, 

(4.11) 

•  W(,  1-  Vn_  i  •  tVj  +  ...  +  v0  •  tv„  =  un. 

Using  the  values  (4.10)  as  initial  approximations,  the  method  of  iterative  residual  correction 
can  then  applied  to  the  system  (4.11)  to  obtain  values  of  tvo,  ti»i, . .  •  ,tvn  €  S  of  maximal 
accuracy  [2j ,  [3j,  [l2j.  Thus,  the  calculation  of  Taylor  coefficients  of  a  function  can  be 
performed  at  the  same  level  of  roundoff  error  as  the  calculation  of  the  value  of  the  function. 
Furthermore,  since  the  Taylor  coefficients  of  the  results  of  arithmetic  operations  are  rational 
functions  of  the  operand  values,  it  is  also  possible  to  extend  the  method  given  by  Bohm 
[2] ,  [3]  for  the  evaluation  of  rational  functions  to  maximal  accuracy  to  the  evaluation  of 
Taylor  coefficients  of  rational  functions  to  maximal  accuracy. 
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